/*
 * @(#)ScaledConjugateGradient.java        1.0 2000/06/26
 *
 * This file is part of Xfuzzy 3.0, a design environment for fuzzy logic
 * based systems.
 *
 * (c) 2000 IMSE-CNM. The authors may be contacted by the email address:
 *                    xfuzzy-team@imse.cnm.es
 *
 * Xfuzzy is free software; you can redistribute it and/or modify it
 * under the terms of the GNU General Public License as published by
 * the Free Software Foundation.
 *
 * Xfuzzy is distributed in the hope that it will be useful, but WITHOUT
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
 * FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
 * for more details.
 */

package xfuzzy.xfsl.model.algorithm;

import xfuzzy.xfsl.model.*;
import xfuzzy.lang.*;
import xfuzzy.xfds.XfdsDataSet;

/**
 * Algoritmo de gradiente conjugado escalado de Moller.
 * 
 * @author Francisco Jos� Moreno Velo
 *
 * @see "M�ller, M.F., A Scaled Conjugate Gradient Algorithm for Fast 
 * Supervised Learning, Neural Networks, Vol. 6, pp. 525-533, 1993."
 */
public class ScaledConjugateGradient extends XfslGradientBasedAlgorithm {

	//----------------------------------------------------------------------------//
	//                            MIEMBROS PRIVADOS                               //
	//----------------------------------------------------------------------------//

	private static final double SIGMA = 1.0e-4;
	
	private static final double LAMBDA = 1.0e-6;

	//----------------------------------------------------------------------------//
	//                            MIEMBROS PRIVADOS                               //
	//----------------------------------------------------------------------------//

	private double lambda;
	
	private double lambdab;
	
	private double[] lp;
	
	private int iter;

	//----------------------------------------------------------------------------//
	//                                CONSTRUCTOR                                 //
	//----------------------------------------------------------------------------//

	/**
	 * Constructor
	 */
	public ScaledConjugateGradient() {
	}

	//----------------------------------------------------------------------------//
	//                             M�TODOS P�BLICOS                               //
	//----------------------------------------------------------------------------//

	/**
	 * Devuelve el c�digo de identificaci�n del algoritmo
	 */
	public int getCode() {
		return SCALED;
	}

	/**
	 * Actualiza los par�metros de configuraci�n del algoritmo
	 */
	public void setParameters(double[] param) throws XflException {
		if(param.length != 0) throw new XflException(26);
	}

	/**
	 * Obtiene los par�metros de configuraci�n del algoritmo
	 */
	public XfslAlgorithmParam[] getParams() {
		return new XfslAlgorithmParam[0];
	}

	/**
	 * Ejecuta una iteraci�n del algoritmo
	 */
	public XfslEvaluation iteration(Specification spec, XfdsDataSet pattern,
			XfslErrorFunction ef) throws XflException {
		XfslEvaluation eval=null;
		XfslEvaluation prev = computeErrorGradient(spec,pattern,ef);
		Parameter[] param = spec.getAdjustable();
		double[] val = new double[param.length];
		double[] g = new double[param.length];
		double[] lg = new double[param.length];
		double[] s = new double[param.length];
		double[] p;
		boolean reset = false;

		if(init) { init=false; iter=0; reset=true; lambda=LAMBDA; lambdab=0;}
		else iter++; 
		if(iter == param.length) { reset = true; iter = 0; }

		for(int i=0; i<param.length; i++) {
			g[i] = -param[i].getDeriv();
			lg[i] = -param[i].getPrevDeriv();
		}

		if(reset) p = g; else p = searchDirection(g,lg,lp);

		double module=0;
		for(int i=0; i<p.length; i++) module += p[i]*p[i];
		double sigma = SIGMA/Math.sqrt(module);
		for(int i=0; i<param.length; i++) val[i] = param[i].value;
		for(int i=0; i<param.length; i++) param[i].setDesp(sigma*p[i]);
		spec.update();
		for(int i=0; i<param.length; i++) param[i].setDeriv(0);
		computeErrorGradient(spec,pattern,ef);
		for(int i=0; i<param.length; i++) s[i]=(param[i].getDeriv() + g[i])/sigma;
		double delta=0;
		for(int i=0; i<p.length; i++) delta += p[i]*s[i];

		boolean success = false;
		while(!success) {
			delta += (lambda - lambdab)*module;
			if(delta<=0) {
				lambdab = 2*(lambda - delta/module);
				delta = -delta + lambda*module;
				lambda = lambdab;
			}
			double mu = 0;
			for(int i=0; i<p.length; i++) mu += p[i]*g[i];
			double alpha = mu/delta;

			for(int i=0; i<param.length; i++) param[i].value = val[i];
			for(int i=0; i<param.length; i++) param[i].setDesp(alpha*p[i]);
			spec.update();
			eval = ef.evaluate(spec,pattern,prev.error);

			double Delta = 2*delta*(prev.error - eval.error)/(mu*mu);
			if(Delta<0) lambdab = lambda;
			else { lambdab = 0; success = true; }
			if(Delta>=0.75) lambda = 0.25*lambda;
			if(Delta<0.25) lambda += delta*(1-Delta)/module;
		}

		for(int i=0; i<param.length; i++) {
			param[i].setPrevDeriv(-g[i]);
			param[i].setDeriv(0);
		}

		lp = p;
		return eval;
	}

	//----------------------------------------------------------------------------//
	//                             M�TODOS PRIVADOS                               //
	//----------------------------------------------------------------------------//

	/**
	 * Calcula la direcci�n de busqueda
	 */
	private double[] searchDirection(double[] g, double[] lg, double[] lp) {
		double mod1=0,mod2=0,mod3=0;
		for(int i=0; i<g.length; i++)
		{ mod1+=g[i]*g[i]; mod2+=g[i]*lg[i]; mod3+=lp[i]*lg[i]; }
		double beta = (mod1 - mod2)/mod3;
		double[] p = new double[lp.length];
		for(int i=0; i<p.length; i++) p[i] = g[i] + beta*lp[i];
		return p;
	}
}

